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Abstract 

This  paper  presents  the  RQ  (Reparameterized  Quaternion)  parameterization  for  three  degree  of  freedom 
(DOF)  rotations  as  an  alternative  to  quaternions  for  the  applications  of  differential  control  and  dynamic 
simulation.  Because  it  requires  just  three  Euclidean  parameters  instead  of  the  four  non-Euclidean  parameters 
of  quaternions,  RQ  is  more  efficient,  less  difficult  to  incorporate  into  large  systems,  and  able  to  make  use  of 
Euclidean  interpolants,  such  as  cubic  splines.  While  RQ  does  possess  singularities  in  its  parameter  space,  we 
show  how  they  can  be  easily  avoided  for  the  intended  applications.  After  discussing  RQ^s  properties  and  limi¬ 
tations  when  applied  to  interpolation  and  spacetime  optimizations,  we  derive  an  extension  to  two  DOF  rota¬ 
tions,  and  show  how  it  can  be  used  to  model  a  ball-and-socket  joint  to  a  very  good  approximation. 
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1  Introduction 

Two  primary  goals  in  computer  animation  are  making  objects  and  characters  move  realistically  and 
controllably.  There  are  many  different  approaches  and  ideologies  for  achieving  these  goals,  but  they  must 
aU  share  a  common  concern:  choosing  a  representation  for  the  underlying  degrees  of  freedom  (DOFs)  of 
the  objects  they  animate.  In  choosing  a  representation,  we  must  consider  three  important  factors: 

•  How  well  does  it  model  the  way  the  object  or  joint  actually  moves? 

•  How  does  it  behave  in  the  intended  application  areas?  This  includes  numerical  properties  and  com¬ 
plexity  of  computing  the  necessary  quantities. 

•  What  level  of  complexity  does  it  impart  to  the  system  we  are  building? 

The  objects  and  characters  being  animated  today  have  many  different  kinds  of  DOFs,  some  of  which 
determine  structural  poses,  and  some  of  which  specify  finer,  more  “organic”  quantities  like  facial  expres¬ 
sions  and  hair.  We  are  concerned  with  the  first  of  these  two  groups,  which  are  known  as  rigid  body 
DOFs,  and  consist  of  translations  and  rotations.  Rigid  body  DOFs  are  the  primary  concern  of  most  com¬ 
mercial  keyframe  animation  systems,  inverse  kinematics  algorithms,  and  ongoing  research  into  dynamic 
simulation  and  optimization. 

Translations  can  be  represented  as  liaear  functions  of  a  set  of  DOFs,  and  are  therefore  straightforward 
to  encode,  understand,  and  use.  Rotations,  however,  are  nonlinear  functions  of  their  DOFs,  which  in¬ 
creases  the  complexity  of  algorithms  needed  to  manipulate  and  control  them.  Further,  the  mathematical 
space  in  which  rotations  are  defined  is  non-Euclidean,  which  complicates  the  numerical  issues  involved  in 
using  derivatives  of  rotational  DOFs,  which  we  must  do  for  inverse  kinematics,  d3mamic  simulation,  and 
optimization. 

The  work  presented  in  this  paper  grew  out  of  a  profound  dissatisfaction  with  available  representations 
of  rotations  for  use  in  our  research,  which  entails  differential  control,  optimization,  and  interpolation  of 
rotations  in  the  context  of  human  character  animation.  We  will  discuss  the  shortcomings  of  the  param- 
eterizations  currently  used  in  graphics  in  section  2.  Then  in  section  3  we  will  present  a  new  parameter¬ 
ization  for  three  DOF  rotations  that,  while  not  without  limitations,  we  believe  represents  a  practical  step 
forward  for  controlling  and  simidating  with  rotations.  In  section  4  we  demonstrate  why  this  is  so,  and 
also  discuss  the  limitations  of  the  parameterization  when  applied  to  interpolation  and  spacetime  optimi¬ 
zation.  In  section  5  we  will  extend  the  parameterization  to  a  two  DOF  rotation  that  is  useful  m  modeling 
ball-and-socket  joints,  and  which  has  few  practical  limitations.  We  also  include  in  the  appendices  formu¬ 
lae  and  C  code  for  computing  and  differentiating  with  the  new  parameterization. 

2  Background 

The  applications  that  are  driving  both  research  and  industry  today  include  inverse  kinematics,  dynamic 
simulation,  and  spacetime  optimization.  Although  these  applications  vary  considerably  in  the  size  and 
complexity  of  the  motion  problems  they  address,  they  all  require  the  same  basic  fimctionality: 

1.  the  ability  to  compute  positions  and  orientations  of  points  on  body  parts  as  fimctions  of  the  parame¬ 
ters 

2.  the  existence  of  and  the  ability  to  compute  derivatives  of  these  positions  and  orientations  with  respect 
to  the  parameters 

3.  the  ability  to  interpolate  smoothly  and  controllably  between  sequences  of  parameter  keyframes 

Regardless  of  the  rotation  parameterization  and  high-level  problem  structure,  rotations  are  generally 
converted  into  3x3  or  4x4  transformation  matrices  as  a  convenient  intermediary  form.  This  gives  us  a 
common  baseline  for  computation  across  various  parameterizations:  we  will  be  interested  in  computing  a 


a)  initial  position  at  0=4=^yt=O  b)  response  to  force  in  +Z  direction 


c)  response  to  force  in  +X  direction  d)  cannot  respond  to  fliis  force 

Figure  1:  Gimbal  lock  in  YXY  Euler  angles.  The  box  can  only  rotate  about  its  center,  and  the  rotation  is  param¬ 
eterized  by  three  Euler  angles  (8,<|),Y),  rotating  first  about  the  Y  axis,  then  the  X  axis,  then  Y  again  (this  particular 
sequence  is  sometimes  used  by  physicists).  When  all  three  angles  are  zero,  the  first  and  last  rotation  axes  are  iden¬ 
tical,  and  it  is  impossible  to  rotate  the  box  about  the  Z  axis,  no  matter  how  great  a  force  we  apply. 


rotation  matrix  and  the  partial  derivatives  of  that  rotation  matrix  with  respect  to  its  parameters.  There¬ 
fore,  if  the  parameterization  possesses  an  n  element  vector  of  DOFs  v  ,  we  must  be  able  to  compute: 

R(v)  and 

dv 

where  R  is  a  4x4  matrix  and  dR/dv  is  a  4x4xm  tensor,  that  is,  an  n  element  vector  of  4x4  matrices,  each 
of  which  is  the  partial  derivative  of  R  with  respect  to  one  of  the  n  parameters  in  v  . 

Interpolating  three  DOF  rotations  is  problematic  because  rotations  are  non-Euchdean,  and  all  of  the 
most  convenient  interpolants,  such  as  the  family  of  cubic  splines,  work  well  only  for  Euclidean  quantities. 
A  best-case  scenario  for  a  parameterization  of  rotations  would  allow  us  to  use  ordinary  cubic  splines  for 
interpolation  of  rotations. 

Now  that  we  know  the  quantities  we  are  interested  in  computing,  we  can  examine  the  parameteriza- 
tions  in  use  today  and  see  why  they  are  imsatisfactory. 

2.1  Euler  Angles 

An  Euler  angle  is  a  DOF  that  represents  a  rotation  about  one  of  the  coordinate  axes.  There  are  three 
distinct  formulae  Rc ,  R;. ,  and  R  for  computing  rotation  matrices,  depending  on  which  coordinate  axis 
the  Euler  angle  represents.  These  formulae  involve  trig  functions  of  the  Euler  angle,  and,  although  they 
are  nonlinear,  their  partial  derivatives  are  easy  to  compute. 

The  Euler  angle  is  really  a  one  DOF  rotation,  and  in  this  capacity  it  functions  well  in  all  respects.  In 
practice,  when  we  need  a  two  or  three  DOF  rotation,  we  replace  R  with  a  sequence  of  single  DOF  Euler 
rotation  matrices.  There  are  several  such  sequences  that  can  be  used  to  represent  a  three  DOF  rotation, 
and  this  in  itself  is  problematic,  since  converting  from  one  to  the  another  is  nontrivial. 

This  concatenation  of  single  DOF  rotations  causes  more  serious  problems  because,  no  matter  what 
sequence  we  choose,  there  will  always  be  at  least  one  configuration  {i.e.  values  for  each  of  the  Euler  an- 


gles)  where  the  first  and  last  axes  of  rotation  in  the  sequence  align,  which  causes  two  of  the  three  partial 
derivatives  to  become  linearly  dependent,  resulting  in  a  loss  of  a  rotational  degree  of  freedom.  This  situa¬ 
tion,  called  gimbal  lock,  means  that  when  in  such  a  configuration,  there  is  a  direction  in  which  we  can 
pull  or  push  with  infinite  force  and  yet  be  unable  to  induce  a  rotation  because  no  combination  of  differen¬ 
tial  change  of  the  three  Euler  angles  will  produce  any  rotation  in  that  direction.  This  situation  is  illus¬ 
trated  in  Figure  1  for  the  sequence  YXY  of  Euler  angles,  which  is  capable  of  representing  any  orientation, 
but  is  gimbal  locked  at  (0,0,0). 

Interpolation  of  multiple  DOF  rotations  is  also  a  problem  for  Euler  angles  because  rotations  confoimd 
axes  (see  [5]  for  more  details),  and  Euler  angles  ignore  this  fact.  When  interpolating  sequences  of  orien¬ 
tations,  the  Euler  angle  model  interpolates  around  each  of  its  three  axes  simultaneously  and  independ¬ 
ently.  This  shows  up,  especially  for  freely  spinning  objects,  as  strange  contortions  and  gyrations. 

In  short,  Euler  angles  are  not  a  good  model  for  multiple  DOF  rotations,  because  they  represent  such 
rotations  as  a  sequence  of  single  DOF  rotations  about  coordinate  axes.  Despite  the  persistent  presence  of 
Euler  angles  in  commercial  ke3dfame  systems,  they  are  imsuitable  for  inverse  kinematics,  d3Tiamic 
simulation,  and  spacetime  optimization  of  three  DOF  rotations. 


2.2  Quaternions 

Euler  himself  showed  that  any  orientation  can  be  achieved  by  a  single  rotation  about  the  proper  axis  (in 
general  not  one  of  the  coordinate  axes).  Quaternions,  as  presented  to  the  graphics  community  by  Shoe- 
make  [5],  can  be  used  as  an  elegant  implementation  of  this  model  that  also  produce  smooth,  natiiral- 
looking  interpolations  of  sequences  of  orientations.  Quaternions  are  4D  objects  that  are  thought  of  as 
having  a  3D-vector  component  and  a  scalar  component.  If  the  quaternion  has  unit  length,  then  the  scalar 
component  can  be  interpreted  as  the  cosine  of  V2  the  angle  of  rotation,  and  the  vector  as  the  axis  of  rota¬ 
tion  scaled  by  sine  of  ^  the  angle  of  rotation,  i.e. 

5  =  [smy^^)ax,  cosO^^)] 

where  9  is  the  angle  of  rotation  and  ax  is  a  unit  length  axis  of  rotation. 

The  formula  for  computing  a  rotation  matrix  R  from  the  unit  quaternion  is  simple  and  inexpensive^,  as 
are  the  partial  derivatives  of  R  with  respect  to  each  of  the  four  quaternion  components.  Fmrther,  the 
partial  derivatives  exist  and  are  linearly  independent  over  the  entire  surface  of  the  imit  4-sphere  on 
which  quaternions  are  defined,  which  means  that  unit  quaternions  are  free  from  gimbal  lock  when  used 
to  control  orientations. 

However,  this  behavior  does  not  come  without  a  cost.  A  quaternion  uses  four  parameters  to  represent  a 
three  DOF  rotation,  which  means  that  there  is  always  a  direction  of  change  (i.e.  partial  derivative)  for  the 
quaternion  that  induces  a  transformation  that  is  not  a  rotation.  This  direction  is,  in  fact,  along  the  4- 
vector  defined  by  the  quaternion,  since  any  change  along  this  direction  would  cause  the  length  of  the 
quaternion  to  become  non-unit,  and  thus  no  longer  a  rotation. 

Several  strategies  have  been  developed  to  deal  with  this  problem.  The  obvious  approach  is  to  add  ex¬ 
plicit  constraints  that  force  quaternions  to  maintain  their  unit  length.  If  the  intended  application  sup¬ 
ports  nonlinear  constraints,  this  can  be  accomplished  with  a  single  scalar  constraint  for  each  quaternion, 
which  forces  the  quaternion  to  maintain  its  imit  length.  While  this  does  work  for  differential  control  and 
simulation,  it  forces  the  systems  of  equations  that  such  applications  solve  to  be  larger  than  necessary, 
since  we  have  an  extra  parameter  and  an  extra  constraint  for  each  rotation. 

Gleicher  [3]  eliminates  the  need  for  the  extra  constraints  by  building  autonormalization  into  the  for¬ 
mula  for  computing  R  from  q .  In  this  scheme,  q  can  take  on  any  value  except  (0,0, 0,0),  and  the 


1  The  rotation  matrix  M  presented  in  [5]  transforms  row  vectors,  so  is  actually  the  transpose  of  the  matrix 
R  that  we  desire  for  transforming  column  vectors. 
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components  of  q  are  replaced  by  (analytically)  normalized  versions  in  the  formulae  for  computing  R  and 
its  four  partial  derivatives.  The  result  is  that  differential  change  in  the  quaternion  can  no  longer  cause  it 
to  cease  representing  a  rotation;  however,  the  quaternion  can  still  change  along  the  direction  of  its  length, 
but  doing  so  now  has  no  tangible  effect  whatsoever.  Thus  for  each  quaternion  there  is  a  DOF  that 
induces  no  change  in  configuration.  This  is  unacceptable  for  applications  such  as  dynamic  simulation 
(and  some  control  algorithms)  that  invert  a  quantity  known  as  the  mass  matrix^  because  the  ineffectual 
DOF  causes  the  mass  matrix  to  become  singular. 

Baraff  and  Witkin  [1]  take  a  different  approach,  discarding  the  quaternion  parameters  as  the  basis  for 
differentially  controlling  a  quaternion  rotation.  They  reason  as  follows:  given  a  current  orientation  repre¬ 
sented  by  a  quaternion  q ,  the  orientation  can  differentially  change  by  acquiring  an  angular  velocity  m  , 
which  is  a  3-vector  whose  direction  gives  an  axis  of  rotation  and  whose  magnitude  is  the  rate  of  rotation 
in  radians  /  time  xmit  about  that  axis.  The  underlying  representation  for  the  rotation  is  still  a 
quaternion,  and  the  formula  for  computing  R  is  the  same,  but  we  no  longer  compute  the  four  partial 
derivatives  of  R  wrt.  q ,  Instead,  we  compute  three  matrices  that  specify  how  the  orientation  will  change 
given  a  differential  angular  velocity  with  components  along  each  of  the  coordinate  axes.  This  formulation 
is  free  from  gimbal  lock,  and  is  also  quite  efficient,  since  there  are  no  extra  constraints  required,  and  the 
number  of  "DOFs”  that  appear  in  the  equations  of  motion  is  three  for  each  rotation. 

There  is,  however,  a  heavy  cost  in  code  complexity,  since  the  number  of  parameters  defining  the 
rotation  (foxir)  is  different  from  the  number  of  DOFs  we  use  to  control  it  (three).  Baraff  provides  formulae 
for  computing  the  time  derivative  of  a  quaternion  from  a  given  angular  velocity,  which  must  be  performed 
before  integrating  the  differential  equations  that  often  arise  in  the  intended  applications.  To  make  the 
most  efficient  use  of  this  technique,  the  disparity  between  nmnber  of  parameters  and  number  of  DOFs 
must  be  propagated  through  large  code  segments.  Furthermore,  it  cannot  be  made  to  work  for 
optimization,  so  its  use  is  limited  to  differential  control  and  simulation. 

Recall  now  that  in  addition  to  derivatives,  we  are  also  interested  in  interpolation  of  orientations.  The 
surface  of  the  unit  4-sphere  really  is  the  correct  place  to  interpolate  rotations,  because,  with  antipodes 
identified,  it  possesses  the  same  metric  and  group  structure  as  SO(3) ,  the  rotation  group.  Shoemake  [5] 
constructs  spherical  Bezier  curves  on  the  surface  of  the  unit  4-sphere,  and  Barr  et  al.  [2]  use  optimization 
techniques  to  solve  for  a  series  of  interpolating  quaternions  that  minimizes  tangential  acceleration  on  the 
surface  of  the  4-sphere  while  satisfying  angular  velocity  constraints  at  the  key  orientations.  However, 
neither  technique  admits  an  analytical  closed  form  (which  precludes  altogether  their  use  as  basis  func¬ 
tions  for  representing  time-varying  orientations),  and  the  quantities  available  for  controlling  the  shape  of 
the  interpolation  are  difficult  to  visualize. 

In  summary,  quaternions  are  quite  well  behaved  for  most  possible  applications,  but  they  incur  an  over¬ 
head  in  efficiency  and/or  code  complexity.  Furthermore,  the  available  means  of  interpolation  are  not 
easily  controlled. 

3  RQ  Parameterization 

What  we  would  really  like  is  a  parameterization  that  behaves  like  a  quaternion  but  is  embedded  in 
Euclidean  space  such  that  the  number  of  parameters  equals  the  number  of  DOFs.  Such  a  parameteriza¬ 
tion  would  be  simple  to  control  and  allow  us  to  use  ordinary  cubic  splines  to  interpolate  rotations.  This 
goal  is,  of  course,  unrealizable,  as  it  is  a  standard  exercise  in  a  topology  text  to  show  that  SO(3)  cannot  be 
mapped  into  without  singularities,  i.e.  gimbal  lock.  Our  more  realistic  goal  is  to  trade  off  some  of  the 
robustness  of  quaternions  for  a  Euclidean  embedding,  trying  to  minimize  the  loss  with  respect  to  the 
mathematical  quantities  we  need  to  compute. 

We  have  devised  a  parameterization  called  RQ  (Reparameterized  Quaternion)  that  meets  these  criteria. 
RQ  represents  a  rotation  by  a  3-vector  v  ,  whose  direction  gives  the  axis  of  rotation  and  whose  magnitude 
is  the  amoimt  of  rotation  about  that  axis  (in  radians).  There  are  several  known  methods  of  translating 
such  an  axis! angle  parameterization  directly  into  a  rotation  matrix,  such  as  Rodrigues*  formula  and 
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Cayley's  parameterization  [4].  However,  we  use  quaternions  as  an  intermediate,  since  they  are  computa¬ 
tionally  equivalent  to  the  direct  methods  that  possess  the  properties  we  desire,  and  are  easier  to  analyze; 
the  formulae  for  computing  the  individual  quaternion  elements  as  functions  of  v  are  as  follows: 

Let 

=  H 

Then 

^ — V,  cosO^O) 

& 


There  are  two  immediate  observations  to  make  about  this  parameterization.  First,  the  quaternion  will 
always  have  unit  magnitude  regardless  of  the  value  of  v ,  which  means  every  point  in  RQ  parameter 
space  represents  a  rotation  -  we  need  take  no  special  measures  to  insure  the  quaternion  has  unit  length. 
Second,  there  seems  to  be  a  singularity  in  the  parameter  space  at  v  =  [0, 0, 0] ,  since  at  this  point  0  is  zero. 
However,  sm(X6)/6  is  an  indeterminate  form^  which  has  continuous  value  and  derivative  in  the  limit  as  0 
goes  to  zero  (in  fact,  it  is  the  well-known  sine  function).  We  can  compute  the  value  at  the  limit  using 
VHopitaVs  Rule,  which  states: 

Suppose  a(;c)  and  b(;c)  are  two  functions  that  are  continuously  differentiable  in  the  neighborhood  sur¬ 
rounding  (but  not  including)  the  point  d  ,  and  b(:r)  is  nonzero  in  die  neighborhood.  Suppose  fiirther  that 

lim  a(x)  =  0  =  lim  b(x)  If  so,  then: 

x^d  x-^d  x-^d  b(jc)  bX^c) 

Since  the  conditions  are  met,  we  can  use  FHopital’s  Rule  to  compute  the  vector  part  of  ^  in  the  limit  as 
0  goes  to  zero  by  differentiating  with  respect  to  0: 

This  gives  the  expected  result  that  v  =  [0, 0, 0]  maps  to  the  zero  rotation  quaternion  q  =  [0, 0, 0,1] .  Fur¬ 
thermore,  away  from  zero,  the  THopitaFs  Rule  formulae  are  numerically  indistinguishable  finm  the  origi¬ 
nal  formulae  out  to  an  appreciable  fraction  of  a  radian.  Therefore,  as  long  as  we  keep  the  transition  point 
from  using  one  set  to  the  other  reasonably  small  (say,  square  root  of  machine  precision),  q  will  be  a  con¬ 
tinuous  fimction  of  v  to  within  machine  precision. 


3.1  Derivatives 

Since  we  are  actually  reparameterizing  quaternions,  we  can  compute  the  derivatives  of  R  with  respect 
to  its  RQ  parameters  by  applying  the  chain  rule  of  differentiation.  That  is,  we  have  R(^(v))  and  wish  to 
compute  3R/dv  ,  which  we  can  compute  as  the  product: 

dv  dq  dv 

Since  we  already  know  how  to  compute  the  partial  derivatives  dR/dq ,  the  only  new  quantities  we  need 
are  the  twelve  partial  derivatives  of  the  quaternion  with  respect  to  its  RQ  parameters  (here  ), 

which  are  given  in  Appendix  A.  Additionally,  Appendix  C  discusses  the  supplemental  C  source  code  for 
computing  R ,  9R/dv  ,  and  other  quantities  presented  later  in  this  paper. 


3.2  Tradeoffs 
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So  far  the  RQ  parameterization  seems  to  fulfill  all  of  our  requirements,  implementing  the  axis/angle 
rotation  model  in  three  Euclidean  parameters.  Before  we  can  discuss  its  application  to  the  animation 
problems  we  have  talked  about,  we  must  be  clear  about  what  we  have  given  up. 


3. 2 A  Combining  Rotations 

One  nice  feature  of  quaternions  is  that  there  exists  a  multiplication  operator  that  can  be  used  to  com¬ 
bine  rotations.  If  and  are  unit  quaternions,  then  the  combined  rotation  q^  that  is  the  result  of  first 
rotating  by  q^  and  then  by  can  be  calculated  as  qz-qi^qiy  where  'o'  represents  quaternion  multiplica¬ 
tion,  as  defined  in  [5].  RQ  (and  any  pure  axis/angle  parameterization)  possesses  no  analogous  operation. 
We  can  combine  two  RQ  vectors  using  any  combination  of  binary  Euclidean  operators,  but  in  no  case  will 
the  result  (in  general)  be  equivalent  to  combining  the  two  rotations  -  to  do  so  we  would  need  to  convert 
the  RQ  vectors  to  their  corresponding  quaternions,  perform  quaternion  multiplication,  then  convert  back, 
incurring  several  trig  and  one  inverse  trig  fiinctions. 

Fortunately,  this  operation  is  typically  not  needed  in  the  applications  we  have  talked  about.  Rotations 
are  changed  only  by  direct  parameter  manipulation,  or  incrementally,  via  their  derivatives.  When  they 
are  combined,  it  is  usually  in  the  context  of  a  transformation  hierarchy,  where  all  rotations  have  already 
been  converted  into  transformation  matrices. 


3,2.2  Singularities 

For  the  purposes  of  control  and  simulation,  the  principal  advantage  of  quaternions  over  Euler  angles  is 
their  fireedom  from  gimbal  lock.  We  already  know  that  the  RQ  parameterization  must  have  singularities, 
so  if  it  is  to  be  useful,  we  must  locate  all  singularities  and  show  how  they  can  be  avoided  at  a  cost  that  is 
outweighed  by  the  benefits  of  RQ. 

As  previously  mentioned,  gimbal  lock  is  equivalent  to  the  disappearance  or  loss  of  linear  independence 
of  rotational  derivatives.  Therefore  we  can  determine  the  singularities  in  RQ  space  by  analytically  ex¬ 
amining  the  partial  derivatives  9^/dv  .  Doing  so  reveals  that  a  singularity  occurs  whenever  is  an  inte¬ 
ger  multiple  of  2n.  Note  that  this  does  not  apply  to  i9=0  because  the  forms  derived  from  THopital’s  Rule 
are  well  behaved.  This  makes  intuitive  sense,  since  a  rotation  of  2%  about  any  axis  is  equivalent  to  no 
rotation  at  all  -  the  entire  shell  of  points  2k  distant  from  the  origin  (and  4;r,  etc.)  collapses  to  the  identity 
in  SO(3) .  This  singularity  manifests  itself  in  the  following  way:  whenever  I  v  I  =  2;z ,  differentially  chang¬ 
ing  V  along  its  own  axis  induces  continued  rotation  about  v  ,  but  differential  movement  in  any  direction 
perpendicular  to  v  induces  a  velocity  for  v  in  the  tangent  plane  of  the  singularity  shelly  thus  (differen¬ 
tially)  keeping  v  on  the  shell  and  inducing  no  change  in  orientation.  In  other  words,  when  v  hits  one  of 
these  shells,  it  can  differentially  keep  rotating  in  the  direction  it  is  already  going,  but  in  no  other  direc¬ 
tion.  Fortunately,  for  the  applications  of  differential  control  and  physical  simulation,  avoiding  these  sin¬ 
gularities  is  easy,  due  to  two  facts. 

First,  if  we  consider  just  the  RQ  parameter  space  inside  the  first  singularity  shell,  each  unique  orienta¬ 
tion  (except  the  rotation  of  zero  radians)  has  two  possible  parameterizations:  as  a  rotation  of  0  radians 
about  V  ,  and  as  a  rotation  Qf2ii-6  radians  about  (quaternions  have  a  similar  property,  in  that  antipo¬ 
des  q  and  -q  represent  the  same  orientation). 

Second,  both  control  and  simulation  operate  by  moving  forward  through  time  in  small  steps,  and  the 
possible  change  in  DOFs  during  each  step  is  small  (much  less  than  a  radian  for  angular  parameters). 

Using  these  facts  it  is  easy  to  avoid  the  singularity:  at  each  time  step  when  the  rotation  is  queried  for 
its  value  and  derivative,  we  examine  |v|,  and  if  it  is  close  to  2ff,  we  reparameterize  v  by  (l“2;r/[v|)v ,  which 
is  an  equivalent  rotation,  but  with  better  derivatives.  Such  dynamic  reparameterization  could,  in  theory, 
also  be  applied  to  avoiding  gimbal  lock  in  Euler  angles,  but  whereas  here  the  reparameterization  simply 
scales  the  current  parameters,  the  corresponding  operation  on  Euler  angles  involves  switching  the  for- 
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mula  used  to  compute  the  rotation  matrix  and  a  sequence  of  inverse  trig  functions  to  determine  the  new 
parameters  [6]. 

4  Applications 

Now  that  we  have  seen  how  the  RQ  parameterization  works  and  what  its  theoretical  limitations  are,  it 
is  time  to  focus  on  the  reason  for  its  creation:  the  simplification  of  algorithms  that  use  parameterized 
rotations.  Of  the  four  applications  we  have  discussed  in  this  paper  -  control,  simulation,  optimization, 
and  interpolation  -  we  beheve  RQ  to  be  ideally  suited  to  the  tasks  of  differential  control  and  d3uiamic 
simulation.  In  the  following  sections  we  will  explain  why,  and  also  discuss  the  complications  that  arise 
in  appl3dng  RQ  to  interpolation  and  optimization. 


4.1  Differential  Control  and  Dynamic  Simtdation 

One  of  the  motivating  applications  for  this  work  is  differential  control,  which  enables  direct  manipula¬ 
tion  interfaces,  inverse  kinematics,  and  real-time  control  of  robotic  manipulators.  To  control  the  positions 
and  velocities  of  objects  and  end-effectors  of  kinematic  assemblies,  the  only  demands  imposed  by  differen¬ 
tial  control  are  that  the  derivatives  dR/dv  be  continuous  and  free  from  gimbal  lock.  Since  control  is  only 
performed  at  discrete  instants  in  time,  the  simple  dynamic  reparameterization  technique  presented  in 
section  3.2.2  wiU  assure  that  these  demands  are  always  met. 

In  dynamic  simulation  applications,  we  track  not  only  an  object’s  instantaneous  position  and  pose,  but 
also  its  linear  and  angular  velocity.  linear  velocity  is  stored  as  a  3-vector  that  represents  the  Cartesian 
direction  and  magnitude  of  the  velocity.  Angular  velocity  is  also  represented  as  a  3-vector  6),  whose 
meaning  is  nearly  identical  to  that  of  an  RQ  vector,  except  that  its  magnitude  0  represents  the  rate  of 
rotation  about  the  axis  rather  than  absolute  orientation,  as  in  RQ. 

To  update  the  position  and  orientation  correctly  as  the  simulation  moves  forward  in  time,  we  need,  in 
addition  to  the  derivatives  necessary  for  differential  control,  a  formula  for  computing  the  time  derivative 
of  the  orientation  given  the  current  orientation  and  the  instantaneous  angular  velocity.  Baraff  [1]  derives 
such  a  formula  for  a  quaternion  orientation: 

q  =  y^a)oq 

where  o)  is  the  angular  velocity  vector  m  extended  with  a  zero  scalar  component  to  make  a  quaternion. 

Since  RQ  does  not  possess  an  equivalent  operation  to  quaternion  multiplication,  we  cannot  derive  a 
similar  formula  for  v .  However,  since  we  can  convert  a  quaternion  into  an  RQ  vector  inside  the  first 
singularity  shell  at  ,  we  can  form  v  in  terms  of  q  like  so: 

dV 

V  =  V(^),  Therefore  v  = - q 

dq 

Y{jq)  involves  one  arctrig  and  one  trig  function,  but  when  we  take  the  derivative  with  respect  to  q  and 
put  the  entire  formula  in  terms  of  v,  much  simplification  occurs.  Appendix  B  gives  a  formula  for  v  in 
terms  of  v  and  (o  that  contains  roughly  the  same  number  and  kind  of  operations  that  the  quaternion 
multiplication  to  compute  q  would  have  required.  Thus  RQ  can  be  used  in  place  of  quaternions  for  dy¬ 
namic  simulation  with  equivalent  performance  and  reduction  in  code  complexity. 

We  claim  that  RQ  is,  in  fact,  ideally  suited  to  these  tasks  because,  despite  the  small  measure  we  must 
take  to  avoid  the  singularity,  implementation  of  RQ  can  be  much  cleaner  than  quaternions.  Not  only  is 
RQ  modestly  more  efficient,  since  it  requires  only  three  DOFs  rather  than  four  and  no  explicit  con¬ 
straints,  but  it  is  also  truly  modular,  in  that  no  constraints  or  information  need  be  propagated  to  code 
segments  outside  the  core  RQ  functions. 
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RQ  Interpolation  Vs.  Spherical  Bezier 
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Figure  2t  Comparison  of  interpolated  rotations  using  RQ  and  spherical  Beziers.  Five  orientations  were  spaced 
evenly  through  four  units  of  time,  and  the  same  tangent  controls  were  used  in  both  cases.  In  the  above  graph,  the 
curves  give  the  interpolated  x,  y,  and  z  values  resulting  from  Euclidean  Bezier  interpolation  of  the  RQ  pa¬ 
rameters;  the  ‘‘Spherical”  curves  are  the  results  of  spherical  Bezier  interpolation  on  the  unit  4-sphere,  converted 
back  into  RQ  parameters. 


4*2  Interpolation 

The  principal  benefit  of  embedding  a  rotation  in  a  Euclidean  space  is  the  possibility  of  using  Euclidean 
interpolants  on  sequences  of  orientations.  Euler  angles  do  not  make  effective  use  of  Euclidean  interpo- 
lants  because  they  ignore  the  interaction  of  rotation  axes.  Since  RQ  is  based  on  the  axis/angle  model  of 
rotations,  we  can  hope  to  do  better. 

However,  the  singularity  shell  (which  is  common  to  all  implementations  of  the  axis/angle  model)  causes 
complications  because  the  distance  between  similar  orientations  grows  without  bound  as  one  approaches 
the  singularity.  We  can  illustrate  this  as  follows:  let  y  =27t-S  for  some  infinitesimal  value  of  S,  Then  if 
Vj  =['>,?,>]  and  Vj  =  ]  are  two  RQ  orientations,  we  note  that  although  they  are  extremely  close  to 

each  other  in  SO(3),  being  only  25  radians  apart,  interpolation  in  RQ  space  will  rotate  through  nearly  An 
radians  -  probably  not  the  desired  behavior.  Quaternions  avoid  this  wraparound  problem  because  inter¬ 
polation  on  the  surface  of  a  sphere  uses  the  same  distance  metric  as  SO(3),  provided  consecutive  orienta¬ 
tions  are  on  the  same  side  of  the  sphere  (t.e.  their  dot  product  is  positive). 

Away  from  the  singularity  shell,  however,  the  results  of  using  corresponding  Euclidean  RQ  interpolants 
and  non-Euclidean  quaternion  interpolants  are  indistinguishable.  Figure  2  shows  a  comparison  of  the 
interpolation  of  five  key  orientations  using  both  spherical  Beziers  on  quaternion  representations  of  the 
orientations  [5]  and  regular.  Euclidean  Beziers  on  RQ  representations.  When  the  quaternions  resulting 
from  the  spherical  Bezier  interpolation  are  converted  into  RQ  parameters,  the  two  methods  are  observed 
to  correspond  quite  closely.  The  fact  that  we  can  achieve  essentially  identical  results  from  interpolating 
in  Euclidean  RQ  space  has  two  important  consequences.  First,  we  can  produce  natxiral  looking  interpola¬ 
tions  using  any  sufficiently  smooth  Euclidean  interpolating  cuive,  including  the  entire  family  of  cubic 
splines,  all  of  which  have  efficiently  computable  closed  forms.  Second,  these  Euclidean  interpolants  have 
fairly  intuitive  controls  in  RQ  space;  for  instance,  if  we  are  using  a  hermite  cubic  interpolant  and  we  want 
an  object  to  have  greater  spin  about  the  y-axis  as  it  is  passing  through  one  of  the  keys,  we  would  simply 
increase  the  y  component  of  the  tangent  control  vectors  at  that  key. 
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So  how  can  we  achieve  this  performance  in  the  face  of  the  problem  near  the  singularity?  In  the  most 
general  case,  when  orientation  keys  are  allowed  to  be  arbitrarily  far  apart  in  parameter  space,  we  can¬ 
not.  But  if  we  stipulate  that  n  consecutive  keys  be  within  one  complete  revolution  of  each  other  (where  n 
is  the  number  of  position  keys  necessary  to  define  a  spline  segment  -  two  for  hermites  and  beziers,  four 
for  b-splines),  then  we  can  dynamically  reparameterize  the  members  (along  with  tangent  controls,  if  they 
exist)  of  each  individual  segment  as  a  group  to  ensure  that  none  are  near  the  singularity  shell,  and  then 
proceed  with  the  interpolation.  This  fix  prevents  us  fi'om  using  cardinal  and  other,  related  splines  that 
cannot  be  broken  up  into  segments.  Also,  to  make  sure  the  stipulation  is  met,  we  may  need  to  insert 
“hidden”  keys  in  between  the  keys  specified  by  the  animator  if  the  change  in  orientation  between  two 
consecutive  keys  is  more  than  one  revolution.  The  need  for  this  last  measure  is,  at  least,  not  exclusive  to 
RQ,  since  quaternions  are  incapable  of  encoding  a  difference  in  orientation  of  An  or  more  between  con¬ 
secutive  keys. 

Finally,  we  note  that  in  practice,  this  added  complication  is  often  not  required.  In  character  animation, 
actors  and  objects  seldom  deviate  from  the  canonical  upright  position  by  more  than  one  revolution  unless 
they  are  in  fireefaU. 


4*3  Spacetime  Optimization 

In  spacetime  and  other  optimizations  that  operate  over  an  entire  animation  simultaneously,  DOFs  are 
not  simply  angles  at  one  instant  in  time,  but  rather  rotation-valued  functions  of  time.  For  iustance,  the 
function  for  a  single  angle  might  be  a  cubic  spline  defined  over  the  animation  time  interval,  in  which  case 
the  DOFs  are  the  positions  of  the  spline^s  control  points,  and  we  need  to  compute  derivatives  of  rotations 
at  various  points  in  time  (i.e.  along  the  curve)  with  respect  to  these  control  points  (simply  several  further 
applications  of  the  chain  rule  to  our  existing  formulae). 

Since  these  functions  are  essentially  interpolating  curves,  we  might  think  about  dealing  with  the  sin¬ 
gularity  in  the  same  manner  as  for  interpolation.  However,  we  cannot  reparameterize  the  functions 
d3mamically,  because  doing  so  will  change  the  shape  of  the  curves,  potentially  perturbing  the  optimiza¬ 
tion  out  of  the  state-space  well  it  is  currently  traversing.  Therefore,  unless  the  possible  range  of  rotations 
is  less  than  2n  (which  actually  is  true  when  solving  for  motion  displacements  rather  than  motions),  we  do 
not  recommend  the  RQ  parameterization  for  spacetime  optimizations  involving  three  DOF  rotations. 


5  Ball-And-Socket  Joints 

The  RQ  parameterization  can  be  even  more  useful  for  modeling  other  types  of  joints  than  freeform  three 
DOF  rotations.  Let  us  consider  the  type  of  ball-and-socket  joint  found  at  the  shoulders  and  the  hips.  If 
we  were  to  map  out  the  legal  range  of  motion  of  these  joints  in  configuration  space  using  RQ  or  any  other 
parameterization,  we  would  get  a  complex,  non-regular  hyper-surface.  However,  we  can  create  a  joint 
that  approximates  the  actual  motion  and  limits  of  motion  fairly  well  using  the  following  scheme,  illus¬ 
trated  in  Figure  3. 

We  model  the  motion  of,  say,  an  arm,  as  a  single  DOF  twist  about  the  principal  axis  of  the  outstretched 
arm,  followed  by  a  two  DOF  swing  of  the  arm  that  incurs  no  twist  about  the  arm^s  axis.  The  last  is  impor¬ 
tant  because  we  wish  to  place  limits  on  the  allowable  twist  of  the  arm;  this  would  be  very  difficult  if  both 
rotations  could  generate  twist,  but  since  only  the  first  rotation  can  twist,  the  limits  can  be  easily  enforced 
with  linear  inequality  constraints.  As  for  creating  suitable  limits  for  the  swing  rotation,  we  will  discuss 
more  accurate  schemes  later,  but  for  now  we  will  stipulate  a  maximum  rotation  in  any  direction,  so  that  a 
rigid  arm  wiU  be  able  to  sweep  out  a  circular  patch  on  the  surface  of  a  sphere. 
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Figure  3:  Degrees  of  Freedom  for  Ball-And-Socket  Joint  The  ‘^arm”  pictured  above  uses  the  ball-and-socket 
model  proposed  here  at  its  shoulder  joint  In  the  image  on  the  left,  the  twist  DOF  is  being  exercised,  and  the  arm 
spins  about  its  axis.  On  the  right  the  two  swing  DOFs  are  used  to  make  the  hand  swing  out  a  circle,  starting  at  the 
bottom  position  and  finishing  closest  to  us.  Note  tiiiat  in  this  entire  motion  there  is  no  qiin  about  the  arm’s  axis. 


We  claim  that  this  is  a  good  model  of  ball-and-socket  joints  because  the  amoimt  of  twist  in  such  joints  is 
fairly  independent  of  the  swing  angle,  and  the  sweep  limits,  while  never  really  circular,  are  often  rather 
ellipsoidal.  Although  we  are  constructing  a  single  joint  from  a  sequence  of  two  rotations,  there  is  no  pos¬ 
sibility  of  gimbal  lock,  since  the  axes  of  rotation  are,  by  definition,  mutually  orthogonal. 

The  twist  can  be  parameterized  as  either  an  Euler  angle,  if  the  principal  axis  for  the  limb  happens  to  be 
a  coordinate  axis,  or  as  a  single  DOF  RQ  rotation,  whose  derivation  is  strai^tforward  given  the  following 
presentation  of  two  DOF  rotations.  Parameterizing  the  twist-less  swing  rotation  would  be  tricky  using 
Euler  angles,  but  given  an  axis/angle  model  like  RQ  it  is  simple,  since  a  necessary  and  sufficient  condition 
for  a  rotation  that  contains  no  spin  about  a  specified  vector  is  for  the  rotation  axis  to  be  orthogonal  to  the 
vector. 

This  means  that  to  achieve  our  desired  two  DOF  rotation,  all  we  need  is  an  RQ  rotation  vector  that  lies 
in  the  plane  perpendicular  to  the  major  axis  of  the  limb  in  its  canonical  (zero  rotation)  position.  We  coidd 
do  this  with  an  explicit  constraint,  but  we  can  more  elegantly  create  such  a  rotation  by  reparameterizing 
the  RQ  rotation  itself  like  so:  pick  any  two  orthogonal,  unit  vectors  in  the  perpendicular  plane;  these 
vectors,  s  and  t  become  a  2D  basis  for  our  desired  swing  rotation  v ,  an  RQ  vector  that  we  compute  like 
so: 

v  =  <xs+Bt 

where  a  and  fi  are  now  the  two  DOFs,  which  we  can  form  into  a  2D  RQ  vector  that  we  will  call  r .  Now 
since  s  and  t  are  unit  vectors,  the  length  of  r  is  the  magnitude  of  the  swing  rotation.  So  to  place  a  limit 
on  the  swing  we  need  only  place  an  inequality  constraint  on  this  vectors  magnitude: 

|r|  <  limit 

In  fact,  it  is  not  much  more  difficult  to  impose  a  more  accurate,  ellipsoidal  angular  limit.  If  we  know  the 
angular  limits  of  rotation  along  the  major  and  minor  axes  of  the  ellipse,  a  and  b ,  then  we  must  choose  s 
and  t  to  correspond  to  the  major  and  minor  axes  in  the  plane,  and  pose  the  following  inequality  con¬ 
straint  instead  of  the  one  above: 
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Since  the  two  DOF  rotation  is  formed  from  a  reparameterization  of  the  three  DOF  rotation,  all  the  same 
formulae  apply,  but  now  we  are  computing  the  3D  RQ  vector  v  from  the  2D  RQ  vector  r ,  To  compute  the 
derivatives  of  R  with  respect  to  r ,  we  simply  apply  the  chain  rule  once  more: 


3R  _  9R  dv 
dr  dq  dv  dr 


where  ^  =  s  and 


This  two  DOF  rotation  is  suitable  for  all  the  same  applications  as  the  three  DOF,  and  additionally  op¬ 
timization.  Recall  that  the  three  DOF  rotation  could  not  be  used  for  optimization  because  dynamically 
reparameterizing  the  control  points  to  avoid  the  singularity  shells  was  imacceptable.  But  the  first  shell 
occurs  at  an  angular  magnitude  of  2;r,  and  since  the  angular  motion  limits  for  the  t3rpes  of  joints  the  two 
DOF  rotation  models  should  always  be  much  less  than  2jl,  crossing  any  of  the  shells  should  never  be  a 
problem. 


6  Conclusion 

We  have  presented  a  new  parameterization  of  three  and  two  DOF  rotations.  Although  it  is  technically 
not  as  robust  as  the  straight  quaternion  parameterization,  it  performs  well  in  practice.  In  practical  terms 
of  flexibihty,  ease  of  interaction,  and  ease  of  integration  into  large  systems,  we  have  found  RQ  to  outper¬ 
form  quaternions  and  all  other  parameterizations  discussed  here. 

At  the  beginning  of  this  paper,  we  stated  three  evaluation  criteria  by  which  any  parameterization  must 
be  judged.  We  will  conclude  by  measuring  RQ  against  them. 


Accuracy  in  modeling  motion.  RQ  is  a  faithful  implementation  of  the  axis/angle  model  of  rotations.  As 
such  it  is  both  a  complete  representation  of  orientations,  and  easily  extendible  to  various  forms  of  re¬ 
stricted  rotations.  We  have  demonstrated  this  by  deriving  a  two  DOF  variant  that  is  quite  useful  in 
modeling  a  ball-and-socket  joint. 


Behavior.  Since  RQ  possesses  singularities  in  its  parameter  space,  it  will  experience  gimbal  lock  at 
complete  revolutions  about  any  axis.  We  have  shown,  however,  how  this  situation  can  be  easily 
avoided  for  most  potential  applications,  excluding  optimization.  Derivative  computation  using  RQ  en¬ 
tails  chain  rule  evaluation,  but  we  have  provided  code  (discussed  in  Appendix  C)  that  encapsulates 
the  computation  efficiently. 


Code  Complexity.  Using  RQ  has  simplified  our  own  code  considerably  over  various  strategies  of  using 
quaternions.  The  "unclean”  aspects  of  RQ  (two  sets  of  formulae  for  computing  with  it;  the  need  for  dy¬ 
namic  reparameterization)  can  be  safely  modularized  and  do  not  propagate  beyond  the  rotation  object. 
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Appendix  A:  Derivatives 

As  in  section  3,  let  v  be  a  3-vector  that  RQ  parameterizes  a  quaternion  q  =  \_q:,y  qy^  q^y  q^  ]  (where  ^^is 
the  scalar  part),  and  9=  [v| 

Now  also  let  i  range  over  the  three  components  of  q  that  make  up  its  vector  part,  and  let  j  range  over 
the  components  of  v  . 

The  formulae  for  computing  the  partial  derivatives  of  q  with  respect  to  v  are,  in  the  usual  case  where 
0»  0: 


^  , 

dVy 


sm(|^) 
'  9 


When/=  j : 
When  19^  j: 


hi  _  1..J  cos(45)  J  sin(|e)  ,  sinde) 
av,  “  ^ 


e" 


,  cos(|^) 


6‘ 

sinCI^) 
'  9^ 


In  the  limit  as  ►  0  and  we  use  rHopitafs  Rule  to  compute  q ,  we  use  the  following  forms  to  compute 
the  partial  derivatives: 


^  =  -^vjcosi^e) 


When  I  =  j: 
When  i  ^  j\ 


h 

av, 


\  cosine) 
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Appendix  B:  v{v,(o) 

Assume  the  same  definitions  as  Appendix  A. 

In  the  normal  case  where  0»  0  (note  that  cot(T0)  =  cos(i0)/sin(-l-0),  and  we  should  have  already  com¬ 
puted  both  these  functions): 

Let  p  =  vxo}  Y  =  9cot(i9) 

0  \  0 

Then  v=:*^(/fi)-;7v  +  p) 

In  the  limit  as  0  ^  0,  the  formula  simplifies  to:  v  =  cci 


Appendix  C:  Sample  Code 

Sample  C  source  code  for  computing  the  rotation  matrix  R ,  its  partial  derivatives  dR/dv  ,  and  v(v,  co)  for 
the  two  and  three  DOF  versions  of  RQ  rotations  can  be  fotmd  at  http://v^vw.cs.cmu .edu/-spiff^rq.  This 
code  can  be  used  to  start  computing  with  RQ  immediately,  but  as  it  stands,  it  is  not  very  efficient.  We 
have  purposely  omitted  the  caching  necessary  to  make  the  implementation  efficient,  because  the  form  of 
cache  storage  will  depend  on  the  modularization  strategy. 

In  particular,  the  values  theta,  cosp ,  sinp,  and  the  quaternion  q  should  be  cached  and  stored  with 
V  when  the  rotation  matrix  is  calculated  from  v,  to  be  used  later  when  calculating  derivatives.  Also,  as 
noted  in  the  code,  when  computing  the  derivatives  of  the  two  DOF  rotation,  the  intermediate  derivatives 
of  q  with  respect  to  a  three  DOF  RQ  rotation  should  be  computed  only  once  and  cached. 
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